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I. INTRODUCTION 



It has been a challenge to understand the relatively high abundances of baryons and anti- 
baryons produced in high energy gold-gold collisions at the Brookhaven National Laboratory 
RHIC (Relativistic Heavy Ion Collider). The center of mass collision energies have ranged 
between 56 and 200 GeV per nucleon pair. This energy is so high that baryon/anti-baryon 
pairs can readily be created. To very good accuracy the relative abundances, including 
multiply strange hyperons, are consistent with them being in chemical equilibrium at a 
temperature of 170 ± 10 MeV and a baryon chemical potential on the order of tens of 
MeV [1,2]. Furthermore, the ratio of baryons to mesons grows with increasing transverse 
momentum, reaching one at px ~ 2 GeV/c [3]. A frequently given explanation of the first 
result is that the hadrons are created in chemical equilibrium at some temperature and 
chemical potential, possibly at the end of a quark-gluon to hadron phase transition. In the 
absence of detailed information on the dynamics it is quite reasonable to postulate that phase 
space is filled randomly. For a big system, this is equivalent mathematically to momentary 
or instantaneous thermal and chemical equilibrium. A possible explanation of the second 
result is that since the transverse momentum is proportional to mass, collective fluid-like 
flow velocity of the expanding matter would boost heavier particles to higher transverse 
momentum. A nice overview of these results can be obtained from the Proceedings of 
Quark Matter 2002 [4]. 

The results on baryon production at RHIC are exciting and stimulating. A parametriza- 
tion of them in terms of temperature T, baryon chemical potential [i arising from the initial 
baryon number of the gold nuclei, and transverse flow velocity vt is a first step. The next 
step is to understand them in terms of microscopic, dynamical processes. In this paper we 
set up a network of rate equations to describe the production of baryon/anti-baryon pairs 
(such as pA and EA) including all members of the octet and decuplet. We also include 
transformations among them (such as A + <-> pn°). The rate of production of octet-octet 
pairs was derived by Kapusta and Shovkovy [5]. To compare with experimental data it is 
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important to include the production of octet-decuplet and decuplet-decuplet pairs too. Un- 
fortunately there is insufficient experimental information to determine these as accurately. 
Therefore, we use a simple form of SU(3) flavor symmetry to estimate these rates. Then 
we solve the network of rate equations in a 2+1 dimensional hydrodynamical model of the 
expanding high energy density matter in heavy ion collisions which has been tuned to repro- 
duce the measured pion multiplicity [3]. Since this is a set of coupled differential equations 
initial conditions are required. We begin to solve the rate equations in the mixed phase of 
a first order quark-gluon to hadron phase transition with three values of T c , namely, 165, 
180, and 200 MeV. In one scenario we assume that all baryons are in chemical equilibrium 
initially, and in the other we assume no baryons initially. Other scenarios are certainly pos- 
sible. We plot the resulting abundances as a function of the local (freeze-out) temperature 
Tf and compare to the measured numbers of p, A, S~, and Q in central gold-gold collisions 
at mid-rapidity at a beam energy of 65 GeV per nucleon at RHIC. We also study what 
happens when the octet-decuplet and decuplet-decuplet rates are increased or decreased by 
a factor of two. Comparison to the transverse momentum distributions of these baryons will 
be done in a later paper. 



II. NETWORK OF RATE EQUATIONS 

We need to specify the dynamical variables, measure or calculate the microscopic rates, 
and then solve the resulting network of rate equations. We shall do each of these in turn. 

Various members of the baryon (or anti-baryon) octet and decuplet have been or will 
be measured in central heavy ion collisions at the CERN Super Proton Synchrotron (SPS) 
and at RHIC. Therefore the variables will be all members of those multiplets. There are 
many mesons that are much lighter than the baryons, such as the pion and the p and uj 
vector mesons to name just a few. The baryons are also much heavier than the three lightest 
species of quarks when the latter are given their current quark masses. Therefore it ought to 
be a good approximation to consider the mesons and/or the quarks and gluons as providing 
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a thermal bath in which baryon/anti-baryon pairs can be created or destroyed. Since kinetic 
equilibrium is usually reached much quicker than chemical equilibrium it is reasonable to take 
the kinetic energy distribution of the baryons and anti-baryons as approximately thermal. 
We assume in this paper that only the absolute number of each species deviates from local 
chemical equilibrium. These approximations can be relaxed but, as they are, the calculations 
become increasingly complex. Ultimately one would reach the point where microscopic 
cascade calculations would be necessary, and then one encounters the problem of multi- 
particle initial states in the quasi-localized interactions. That problem is largely avoided in 
the present approximation because of the thermal averaging of the microscopic rates (see 
below) . 

Under the above conditions, the rate equation for the spatial density of a baryon species 

b is 
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The first term on the right involves the rate for producing the specified baryon/anti-baryon 
pair by strong interaction currents, R{bb') = the number of such pairs produced per unit 
volume per unit time, and is the driving term in the rate equation. Of course, the same pair 
can annihilate each other, and this is related to the production rate by detailed balance. The 
factor in square brackets enforces detailed balance. The n^? ml is the equilibrium density of 
the species b' as represented by the temperature and chemical potentials at the time t. The 
second term is a dilution term; the density will decrease in inverse proportion to the volume 
for an expanding system. The third term arises if there exists a baryon species that can 
decay via the strong interactions into the baryon species of interest. The quantity in square 
brackets following it allows for the inverse reaction and satisfies detailed balance. The last 
term arises if the baryon species of interest can decay into another baryon via the strong 
interactions. There is one first order in time, nonlinear, rate equation for each species of 
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baryon. This makes a coupled network of differential equations. 

The decay rate T(b' — > b + X) is just the inverse lifetime for the specified decay, namely, 
1/t(6' — > 6+X). These are taken from the Particle Data Tables [6]. We include the following 
decays: A -> iV + vr, S* -> A + vr, S* -> E + vr, 2* -> S + vr. 

The rate i? is the crucial ingredient in this network of rate equations. It was derived in 
quite some detail by Kapusta and Shovkovy [5] for the case of baryons and anti-baryons in the 
lowest octet. The basis for that calculation involved an effective current-current interaction 
between the baryons and strong interaction currents such as the vector and axial-vector. 
A combination of hadron phenomenology, experimental data on cross sections, and SU(3) 
flavor symmetry allows the rates to be expressed in terms of the spectral densities of the 
strong interaction currents. Since the y/s is so large, with a minimum value determined by 
threshold production of the baryon/anti-baryon pair, these spectral densities can be taken 
from perturbative QCD. For the production of baryons in the octet, the rate was expressed 
as [5] 

RsAbih) = c+ihb^n+ihh) + C_(M 2 )7MM 2 ) . (2) 

Here 

TZ^hh) = 9{ \t a :l 7 l ) 4 T8 z*z 2 2 {Az 1 K 1 (z 1 )K 2 (z 2 ) + Az 2 K x {z 2 )K 2 {z,) 
4 ( 27r ) it 

±( Zl ± z 2 fK x {z x )K x {z 2 ) + [16 + ( Zl ± z 2 ) 2 ]K 2 ( Zl )K 2 (z 2 )} Fl NN (s) (3) 

where Zi — m,i/T,T is the temperature, and -Fann( s ) is an annihilation form factor evaluated 
at the average value s = (mi + m 2 ) 2 + 3(mi + m 2 )T. A simple parametrization of this form 
factor which is consistent with nucleon/anti-nucleon annihilation data was found to be 

-FaNn(s) = r -, ; \2M \2 ' ( 4 ) 

2.21 + [s — (mi + m 2 ) z \/K z 

with A = 1.63 GeV. According to the latest data analysis [6] the strong interaction coupling 
is a s (m 2 ) = 0.35 ± 0.03. The ± correspond to vector /axial- vector contributions to the rate. 
The numerical coefficients C±(bJ) 2 ) are given in [5]; they are all of order one. 



There is unlikely ever to be sufficient experimental information from hadron-hadron 
scattering to pin down the production rates for baryons in the decuplet since they are all 
unstable with lifetimes less than 1CT 10 seconds. Therefore, we have estimated the rates and 
parametrize them as follows. 



We take the coefficients of the vector and axial-vector contributions to be equal which means 
that 



The coefficients C 8) i (6i6 2 ) and Ci ,io(6i6 2 ) are displayed in Tables I and II, respectively. 
We have constructed them as follows. First, an entry in the table is zero if the quantum 
numbers of the pair do not match the quantum numbers of the vector or axial- vector meson 
nonets. Second, exact isospin symmetry within SU(2) multiplets is used. Third, the sum 
of the coefficients of each column of Table I equals 4, which is the order of magnitude of 
the corresponding octet-octet coefficients C 8j8 ; the sum of the coefficients of each row Table 
I is either 4|, 5, or 5|. Flavor SU(3) symmetry cannot give more information than this. 
Similarly the sum of the coefficients of each column (and row) of Table II equals 6. In our 
numerical work we shall study the consequences of increasing or decreasing these coefficients 
overall by a factor of 2. 

There will always be more baryons than anti-baryons present in the final state of a heavy 
ion collision because the colliding nuclei have a net baryon number. However, this asymmetry 
is reduced as the beam energy increases because produced baryons always come with an anti- 
baryon due to baryon number conservation. The pair production is an increasing function 
of beam energy and so the initial number of baryons from the nuclei becomes a smaller and 
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smaller fraction of the total. At the highest energies at RHIC the baryon chemical potential 
is usually estimated to be on the order of 10 MeV when the temperature is 150 MeV. For 
simplicity of calculation and presentation we shall take the baryon chemical potential to be 
exactly zero in this paper. To remove the difference between the experimental data and our 
calculation due to this approximation, we compare our results to the observed average yield 
of a baryon and antibaryon, for example (p + p)/2, not to the observed baryon yield as such. 

III. COMPARISON TO DATA 

Baryon and anti-baryon production cannot begin until the local energy density is low 
enough, that is, when the matter is in the hadronic phase as opposed to the quark-gluon 
phase. This will occur during the expansion stage of a high energy nucleus-nucleus collision. 
This stage is frequently modeled with hydrodynamics [7]. We shall do so in this paper 
too. We will use a 2+1 dimensional description of the final stage expansion that takes into 
account transverse expansion. The details of these calculations have been given many times 
before and so we just refer the reader to those papers [8]. Briefly, the expansion begins in 
the quark-gluon plasma phase with an energy density adjusted to reproduce the measured 
pion rapidity density at mid-rapidity. The phase transition is first order; we have chosen T c 
to be 165, 180, and 200 MeV and test the sensitivity of baryon production to it. 

The network of rate equations is solved within each co-moving cell in coordinate space. 
Different cells evolve somewhat differently in space and time, and this makes for a very 
computationally intense task. At each point in time, as measured by an observer at rest 
in the center-of- momentum frame of the central gold-gold collision, each cell has its own 
temperature. This local temperature enters in the production rates and in the local chemical 
equilibrium densities. In order to concisely display the results of our calculations we plot 
the solutions to the network of rate equations as a function of the local temperature and not 
as a function of the local time. Because each fluid cell has its own temperature at any given 
time, we must integrate over a constant temperature hypersurface to get the baryon rapidity 
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density at this temperature. Thus the value at a fixed temperature does not correspond to 
any particular time. This choice is rather natural since baryon production begins at a system- 
wide value of the temperature equal to T c and ends at a supposed freeze-out temperature 
Tf where the hadrons lose local thermal equilibrium and begin their free-streaming stage. 

In figure 1 we display the results of numerical solution of the rate equation network for p, 
A, S~ and fl. Plotted is the ratio of the calculated density to the equilibrium density at that 
particular temperature. No weak decays of unstable baryons have been allowed for at this 
point. The rate equations require an initial condition and we have chosen two: all baryons in 
chemical equilibrium at T c , and no baryons present at all at T c . The former are represented 
in the figure by the upper set of thick curves while the latter is represented by the lower set 
of thin curves. Obviously this calculation has used T c = 180 MeV, but it is representative 
of other choices. When baryons are present initially with their equilibrium abundances, 
they evolve in such a way that they always stay above the equilibrium abundance at each 
temperature below T c . The reason is that the system expands more rapidly than chemical 
reactions can keep up with. In other words, the typical annihilation rates are smaller than 
the expansion rate. When no baryons are present initially, their abundance at first builds up 
rapidly. The reason that there is a finite abundance at T c already is that the cells generally 
remain within the mixed phase at T c for a finite time span, thus allowing a buildup before the 
temperature drops below T c . In the range from 120 to 130 MeV, these baryon species have 
caught up to the equilibrium abundance at that temperature. Thereafter they are above 
the equilibrium values for the same reason as stated above, namely, that annihilation rates 
are generally smaller than the expansion rate. Naturally, the abundances with no baryons 
present initially never catch up with the abundances where baryons were initially created in 
chemical equilibrium. 

In figures 2-5 we show the rapidity density at mid-rapidity as a function of freeze-out 
temperature for p, A, S~ and Q, respectively. The upper set of curves result from having 
baryons in chemical equilibrium initially, while the lower set of curves result from having no 
baryons present initially. The central curve in each set results from using the standard set of 
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octet-decuplet and decuplet-decuplet coefficients shown in the Tables. The upper and lower 
curves in the lower set result from increasing and decreasing these coefficients by a factor 
of 2 respectively (the octet-octet coefficients are left unchanged). In the upper set of curves 
the effect is reversed: the larger rates lead to smaller yields. Also shown in these figures 
are the experimental data from RHIC experiments. Because we use the approximation of 
zero baryon chemical potential, we compare our result to the measured average yield of 
baryon and antibaryon. The (p + p)/2 and (A + A)/2 data come from PHENIX [9] and the 
(S~ + S+)/2 and (ft + ft)/2 data come from STAR [10]. The calculation and the PHENIX 
data are for 5% most central collisions whereas the STAR data are for 10% most central 
collisions. By comparing the hydrodynamically calculated pion yields at different centralities 
we have estimated that this leads to about 10% larger yield in our calculation compared to 
the data. The darker central band in each figure represents statistical errors only while the 
lighter outer band includes systematic errors as well. (For ft the systematic error is smaller 
than the statistical error.) We have performed decays of unstable baryons appropriate to 
the particular measurement. For example, protons coming from the weak decay A — > p-\-n~ 
are not included in figure 2, but protons from the decay S + — > p + n° are. 

What we learn from this set of figures is that the baryon production rates are too small 
in comparison to the expansion rate to reproduce any of the experimental data when there 
are no baryons present initially. There must be a significant abundance of baryons present 
already at the beginning of the hadronic phase. Assuming the baryons to be produced in 
chemical equilibrium and then evolving them during the expansion of the matter provides 
agreement with the data on p, A, and H~ within systematic error bars for 120 < Tf < 
T c . However, the theoretical calculations produce about 50% more ft than observed. It is 
interesting that this magnitude of discrepancy does not occur for the S~, even though it has 
two valence s-quarks. Oftentimes in purely statistical models an s-quark fugacity factor is 
used to fit the data. A typical value might be 0.9, so that the A yield would be multiplied 
by 0.9, the yield by (0.9) 2 , and the ft yield by (0.9) 3 . Any such suppression factor should 
be an outcome of the results of a rate equation calculation like the present one, and cannot 
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be done arbitrarily at the end of the calculation to fit data. On the other hand, if one had 
confidence that there was a first-order phase transition with this numerical value of T c and 
that the expansion was described adequately by hydrodynamics, then one could adjust the 
initial conditions to match the experimental data. 

Next we test the sensitivity of the results to the critical temperature. The results are 
shown in figures 6-9. In these calculations we have used the set of coefficients from Tables I 
and II. The calculated abundances are generally increasing functions of T c , chosen here to be 
165, 180 and 200 MeV. As seen in figures 2-5, the abundances calculated when no baryons 
are present initially are all below the data. When the baryons are in chemical equilibrium 
initially, the data on A and 2~ are bracketed by the curves. The p are consistent, within 
systematic errors, with 180 and 200 MeV. The f2 data is consistent with 165 MeV but 
is below the 180 and 200 MeV curves. However, there is no value of T c which provides 
agreement between the calculated abundance and the observed one for all four species of 
baryons. 

IV. CONCLUSION 

In this paper we have estimated the thermal production rates of baryon/anti-baryon 
pairs when one or both of them are members of the spin-3/2 decuplet. This estimate 
is based on the thermal production rates when both are members of the spin- 1/2 octet 
[5]. The latter rates are quite well-founded in hadronic phenomenology and QCD, whereas 
the rates involving the spin-3/2 baryons are somewhat uncertain as they are unstable and 
beams cannot be made to measure cross sections. We used exact isospin symmetry within 
multiplets; flavor SU(3) was used too, but it is insufficient to pin down all the rate coefficients. 
We made numerical estimates based on octet-octet production and ultimately tested the 
sensitivity of the results to changes of a factor 2 in the coefficients. 

We used 2+1 dimensional hydrodynamics, which includes transverse expansion, to model 
the expansion stage of a central gold-gold collision at RHIC. This provided the thermal bath 
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within which the network of rate equations was solved. The network of equations requires 
initial conditions. We chose two interesting limits to study in detail: no baryons present 
initially and baryons initially in chemical equilibrium. We also studied the sensitivity of the 
results to the critical temperature of an assumed first order phase transition from quarks 
and gluons to hadrons. 

When comparison is made to RHIC data on p, A, and Q, the most important con- 
clusion is that the calculated production rates are too small to compensate for the rapid 
expansion of the matter as described by hydrodynamics if there are no baryons present at 
the end of the phase transition. There is just not enough time to build up the abundances 
of these baryons starting from nothing. On the other hand, the observed abundances of p, 
A, and H~ can be reproduced when this network of rate equations is solved within a hydro- 
dynamic expansion if baryons are taken to be in chemical equilibrium when each fluid cell 
is converted from quarks and gluons to hadrons. The observed abundance of Q, however, 
is generally less than the calculated one, and when the parameters are chosen to obtain the 
correct Q abundance the other baryons are underproduced. So the fl remains a bit of a 
puzzle. 

The results obtained here are not trivial nor are they obvious. The interplay among 
the initial conditions, the expansion dynamics (in particular the expansion rate), and the 
production rates of octet and decuplet baryons is intricate and subtle. If the baryons are 
produced in chemical equilibrium at T c , their numbers will decrease as the system expands 
and cools by an amount that depends on the annihilation rates compared to the expansion 
rate. If no baryons are present at T c , their numbers will at first increase with time until 
local equilibrium is achieved at some temperature below T c which is generally different for 
each species. After that there will be an overabundance of baryons relative to equilibrium, 
and their numbers will begin to decrease. 

Among the extensions of this project under investigation are different choices of the initial 
abundances and the transverse momentum distributions of the various baryon species. The 
latter would allow us to study the expansion dynamics in finer details than just the overall 
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dN/dy at mid-rapidity. 
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Table I: Relative strengths of coefficients in the 
expressions for the rates. 
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Table II: Relative strengths of coefficients in the 
expressions for the rates. 
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FIG. 1. The ratio of the calculated abundance of the indicated species to the chemical equilib- 
rium value as a function of the local temperature. The upper set of curves start with the baryons 
in equilibrium at T c , the lower set start with no baryons. 
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FIG. 2. The number of protons at central rapidity in 5% most central gold-gold collisions as 
a function of local temperature. The solid curves start with baryons in equilibrium, the lower set 
with no baryons. The central curve in each set uses the coupling coefficients given in the Tables, 
while the other two use values larger and smaller by a factor of two. With the exception of A, 
baryons unstable to strong or weak decays have been decayed. The data are from the PHENIX 
collaboration [9] for the same centrality. The dark band represents statistical and the light band 
systematic errors. 
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FIG. 3. The number of lambdas at central rapidity in 5% most central gold-gold collisions as 
a function of local temperature. The solid curves start with baryons in equilibrium, the lower set 
with no baryons. The central curve in each set uses the coupling coefficients given in the Tables, 
while the other two use values larger and smaller by a factor of two. With the exception of A, 
baryons unstable to strong or weak decays have been decayed. The data are from the PHENIX 
collaboration [9] for the same centrality. The dark band represents statistical and the light band 
systematic errors. 
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FIG. 4. The number of cascades at central rapidity in 5% most central gold-gold collisions as 
a function of local temperature. The solid curves start with baryons in equilibrium, the lower set 
with no baryons. The central curve in each set uses the coupling coefficients given in the Tables, 
while the other two use values larger and smaller by a factor of two. Baryons unstable to strong 
decays have been decayed. The data are from the STAR collaboration [10] for 10% most central 
collisions. The dark band represents statistical and the light band systematic errors. 



1 Q 



1.0 



0.8 



i — i — r 



i — i — r 



i — i — r 



i — i — r 



Q 



t3 



X5 



0.6 



0.4 



Q Q r= l — I — h ^ ~t~ 1 — I — I — I — \ — t- 4= 3= F= l = 

100 120 140 160 180 

T f (MeV) 

FIG. 5. The number of omegas at central rapidity in 5% most central gold-gold collisions as 
a function of local temperature. The solid curves start with baryons in equilibrium, the lower set 
with no baryons. The central curve in each set uses the coupling coefficients given in the Tables, 
while the other two use values larger and smaller by a factor of two. Baryons unstable to strong 
decays have been decayed. The data are from the STAR collaboration [10] for 10% most central 
collisions. The dark band represents statistical errors only; the systematic errors are smaller than 
the statistical ones. 
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FIG. 6. Same as figure 2 except that the calculations are done for three critical temperatures 
of 165, 180, and 200 MeV. The coefficients from the Tables are used. 
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FIG. 7. Same as figure 3 except that the calculations are done for three critical temperatures 
of 165, 180, and 200 MeV. The coefficients from the Tables are used. 
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FIG. 8. Same as figure 4 except that the calculations are done for three critical temperatures 
of 165, 180, and 200 MeV. The coefficients from the Tables are used. 
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FIG. 9. Same as figure 5 except that the calculations are done for three critical temperatures 
of 165, 180, and 200 MeV. The coefficients from the Tables are used. 
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